# install.packages("devtools")
# devtools::install_github("MathiasHarrer/dmetar")


install.packages("dmetar")
install.packages("meta")
install.packages("metafor")
install.packages("corrplot")
install.packages("Hmisc")
library("dmetar")
library("meta")
library("metafor")
library("corrplot")
library("Hmisc")


GLPT1DMMETA_r= read.rm5(
  "GLP-1 for Type 1 Diabetes- metareg.rm5",
  sep = ",",
  quote = "\"",
  "GLP-1 for Type 1 Diabetes- metareg",
  numbers.in.labels = TRUE,
  debug = 0
)

c1o1 <- metacr(GLPT1DMMETA_r,comp.no=1, outcome.no=1)
year1 <- c(2016, 2016, 2016, 2019, 2018, 2016, 2019, 2020, 2018, 2015, 2020, 2012, 2013, 2020, 2020, 2016, 2016, 2016, 2016, 2016, 2016, 2020, 2014, 2021, 2022)
#Location:   MN = 0, NA, Europe, Asia = 1
locna1 <- c(0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 1, 1, 1, 0, 0, 0, 0, 1, 0, 0)
loceu1 <- c(0, 0, 0, 1, 0, 1, 1, 1, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 1)
locas1 <- c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0)

week1 <- c(26, 26, 26, 26, 52, 24, 52, 26, 24, 12, 26, 52, 52, 24, 26, 12, 12, 12, 26, 26, 26, 52, 52, 54, 12)


#Regimen 0= lira 1= exena, albi
exena1 <- c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0)
albi1 <- c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0)
#Liradose 0= exena, albi  0.6, 1.2, 1.8= lira
liradose1 <- c(0.6, 1.2, 1.8, 1.8, 1.8, 1.8, 1.8, 1.8, 1.8, 1.2, 1.8, 0.9, 0, 0, 0, 0.6, 1.2, 1.8, 0.6, 1.2, 1.8, 0, 0, 1.8, 1.2)

#biases:   1= low risk, 2= medium, 3= high risk
rand1 <- c(1, 1, 1, 1, 2, 1, 2, 1, 1, 1, 1, 3, 3, 2, 1, 2, 2, 2, 1, 1, 1, 1, 3, 1, 2)
attr1 <- c(2, 2, 2, 3, 2, 2, 2, 1, 1, 3, 3, 2, 1, 2, 3, 3, 3, 3, 2, 2, 2, 1, 3, 3, 1)
blin1 <- c(1, 1, 1, 1, 2, 1, 2, 1, 1, 1, 1, 3, 3, 1, 1, 1, 1, 1, 1, 1, 1, 1, 3, 1, 1)
adhr1 <- c(2, 2, 2, 3, 2, 2, 2, 1, 1, 3, 3, 2, 1, 2, 3, 3, 3, 3, 2, 2, 2, 1, 3, 3, 1)
detec1 <- c(1, 1, 1, 1, 2, 1, 2, 2, 2, 2, 2, 2, 2, 2, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2)
repoth1 <- c(1, 1, 1, 1, 2, 1, 2, 1, 2, 1, 1, 2, 2, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2)


# cpep1 <- c(15.2, 15.2, 15.2, NA, NA, 21.0, 100, 9.1, NA, 0, NA, NA, 100, 41.7, 22, 0, 0, 0, 17.4, 17.4, 17.4, NA, NA, 100, 100)
# pump1 <- c(25.0, 26.0, 25.0, 17.9, NA, 0, NA, 100, 47, NA, 69, NA, NA, 73.4, 0, 80, 80, 80, 24.0, 28.3, 29.5, 0, NA, NA, NA)

age1 <- c(43.3, 42.8, 43, 50.4, 46.7, 48, 27, 46.5, 35.8, 37.8, 46, 48.5, 27.7, 36.1, 50.2, 47.5, 46, 46, 43.5, 43.7, 43.6, 22.5, 37.3, 29.5, 33.6) 
# sex1 <- c(45, 47, 46, 80, NA, 65, NA, 32, 47, 66, 38, NA, 72, 32, 72.5, 45.7, 50, 34.4, 47.5, 48.2, 48.9, 56.5, 46, 65.3, 50)
dur1 <- c(20.9, 20.9, 21.1, 32.4, 22.3, 22.5, 0.1, 20.5, 20.3, 18.9, 18.6, NA, 0.1, 19.6, 21.1, 19, 20, 20, 21.3, 21.6, 21.6, 0.2, 20.5, 0.25, 3.4)
ins1 <- c(60.3, 60.3, 60.1, NA, NA, 59.5, 22.5, 51, 70.9, 59.6, 57.5, 39.2, NA, 49.5, 60.4, 49.5, 58.7, 47.1, 58.2, 58.5, 58.9, NA, 42.7, NA, 17.2)
bmi1 <- c(28.9, 28.9, 28.9, 29, 28.9, 30.1, 24, 29.5, 30.5, 23.4, 31.6, 22.8, 21.5, 29.3, 28.3, 27, 30.5, 28, 29.7, 29.6, 29.7, 22.3, 26.1, 24.1, 22.9)
a1c1 <- c(8.08, 8.08, 8.08, 8.2, 7.82, 8.7, 8.4, 8.15, 7.4, 8.75, 7.89, 7.65, 9.7, 7.6, 8.25, 7.57, 7.57, 7.57, 8.16, 8.16, 8.16, 7.29, 7, 7.25, 6.8)






year1 <- scale(year1, scale = FALSE)
locna1 <- scale(locna1, scale = FALSE)
loceu1 <- scale(loceu1, scale = FALSE)
locas1 <- scale(locas1, scale = FALSE)
week1 <- scale(week1, scale = FALSE)
exena1 <- scale(exena1, scale = FALSE)
albi1 <- scale(albi1, scale = FALSE)
liradose1 <- scale(liradose1, scale = FALSE)
rand1 <- scale(rand1, scale = FALSE) 
attr1 <- scale(attr1, scale = FALSE) 
blin1 <- scale(blin1, scale = FALSE)
adhr1 <- scale(adhr1, scale = FALSE)
detec1 <- scale(detec1, scale = FALSE)
repoth1 <- scale(repoth1, scale = FALSE)
age1 <- scale(age1, scale = FALSE)
dur1 <- scale(dur1, scale = FALSE)
ins1 <- scale(ins1, scale = FALSE)
bmi1 <- scale(bmi1, scale = FALSE)
a1c1 <- scale(a1c1, scale = FALSE)




factors1 <- data.frame(Publication_Year= year1,
                       Continent_NA_MN= locna1,
                       Continent_EU_MN= loceu1,
                       Continent_AS_MN= locas1,
                       Study_Duration= week1,
                       Med_Exena_Lira= exena1,
                       Med_Albi_Lira= albi1,
                       Liraglutide_Dose= liradose1)
corm1 <- rcorr(as.matrix(factors1))
# corm1$r[corm1$r < -2] <- NaN
corm1
corrplot(corm1$r, method = 'number')


factors1 <- data.frame(Randomization_Bias= rand1,
                       Attrition_Bias= attr1,
                       Performance_Bias= blin1,
                       Nonadherence_Bias= adhr1,
                       Detection_Bias= detec1,
                       Reporting_Other_Bias= repoth1)
corm1 <- rcorr(as.matrix(factors1))
corm1
corrplot(corm1$r,  method = 'number')


factors1 <- data.frame(Randomization_Bias= rand1,
                       Attrition_Bias= attr1,
                       Detection_Bias= detec1,
                       Reporting_Other_Bias= repoth1)
corm1 <- rcorr(as.matrix(factors1))
corm1
corrplot(corm1$r,  method = 'number')



factors1 <- data.frame(Age= age1,
                       Male_Sex= sex1,
                       DM_duration= dur1,
                       Baseline_A1c= a1c1, 
                       Baseline_BMI= bmi1,
                       Daily_Insulin= ins1)
corm1 <- rcorr(as.matrix(factors1))
corm1
corrplot(corm1$r,  method = 'number')


factors1 <- data.frame(DM_duration= dur1,
                       Baseline_A1c= a1c1)
corm1 <- rcorr(as.matrix(factors1))
corrplot(corm1$r,  method = 'number')




factors1 <- data.frame(Publication_Year= year1,
                       Continent_NA_MN= locna1,
                       Continent_EU_MN= loceu1,
                       Continent_AS_MN= locas1,
                       Study_Duration= week1,
                       Med_Exena= exena1,
                       Liraglutide_Dose= liradose1,
                       Randomization_Bias= rand1,
                       Attrition_Bias= attr1,
                       Detection_Bias= detec1,
                       Reporting_Other_Bias= repoth1,
                       DM_Duration= dur1,
                       Baseline_A1c= a1c1)
corm1 <- rcorr(as.matrix(factors1))
# corm1$r[corm1$r < -2] <- NaN
corm1
corrplot(corm1$r, method = 'number')



c1o1d= data.frame(TE= c(c1o1$TE), seTE= c(c1o1$seTE), 
                  Publication_Year= year1,
                  Continent_NA_MN= locna1,
                  Continent_EU_MN= loceu1,
                  Continent_AS_MN= locas1,
                  Study_Duration= week1,
                  Med_Exena_Lira= exena1,
                  Liraglutide_Dose= liradose1,
                  Randomization_Bias= rand1,
                  Attrition_Bias= attr1,
                  Detection_Bias= detec1,
                  Reporting_Other_Bias= repoth1,
                  DM_Duration= dur1,
                  Baseline_A1c= a1c1)

c1o1.mi <- multimodel.inference(TE= "TE", seTE = "seTE", data = c1o1d,
                     predictors = c("Publication_Year", 
                                    "Continent_NA_MN", 
                                    "Continent_EU_MN", 
                                    "Continent_AS_MN", 
                                    "Study_Duration",
                                    "Med_Exena_Lira", 
                                    "Liraglutide_Dose",
                                    "Randomization_Bias",
                                    "Attrition_Bias",
                                    "Detection_Bias",
                                    "Reporting_Other_Bias",
                                    "DM_Duration", "Baseline_A1c"),
                     interaction=FALSE)
c1o1.mi
#
c1o1.mreg <- rma(yi = TE,
                sei = seTE,
                data = c1o1,
                method = "ML",
                mods = ~sex1,
                test = "knha")
c1o1.mreg
c1o1.mreg1 <- rma(yi = TE,
                 sei = seTE,
                 data = c1o1,
                 method = "ML",
                 mods = ~loceu1+liradose1+sex1,
                 test = "knha")
c1o1.mreg1


anova(c1o1.mreg, c1o1.mreg1)



c1o2 <- metacr(GLPT1DMMETA_r,comp.no=1, outcome.no=2)
year12 <- c(2016, 2016, 2016, 2019, 2018, 2016, 2020, 2018, 2015, 2020, 2013, 2020, 2020, 2016, 2016, 2016, 2016, 2016, 2016, 2020, 2014, 2021, 2022)
#Location:   MN = 0, NA, Europe, Asia = 1
locna12 <- c(0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 1, 0, 1, 1, 1, 0, 0, 0, 0, 1, 0, 0)
loceu12 <- c(0, 0, 0, 1, 0, 1, 1, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 1)
locas12 <- c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0)

week12 <- c(26, 26, 26, 26, 52, 24, 26, 24, 12, 26, 52, 24, 26, 12, 12, 12, 26, 26, 26, 52, 52, 54, 12)


#Regimen 0= lira 1= exena, albi
exena12 <- c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0)
albi12 <- c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0)
#Liradose 0= exena, albi  0.6, 1.2, 1.8= lira
liradose12 <- c(0.6, 1.2, 1.8, 1.8, 1.8, 1.8, 1.8, 1.8, 1.2, 1.8, 0, 0, 0, 0.6, 1.2, 1.8, 0.6, 1.2, 1.8, 0, 0, 1.8, 1.2)


#biases:   1= low risk, 2= medium, 3= high risk
rand12 <- c(1, 1, 1, 1, 2, 1, 1, 1, 1, 1, 3, 2, 1, 2, 2, 2, 1, 1, 1, 1, 3, 1, 2)
attr12 <- c(2, 2, 2, 3, 2, 2, 1, 1, 3, 3, 1, 2, 3, 3, 3, 3, 2, 2, 2, 1, 3, 3, 1)
detec12 <- c(1, 1, 1, 1, 2, 1, 2, 2, 2, 2, 2, 2, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2)
repoth12 <- c(1, 1, 1, 1, 2, 1, 1, 2, 1, 1, 2, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2)


dur12 <- c(20.9, 20.9, 21.1, 32.4, 22.3, 22.5, 20.5, 20.3, 18.9, 18.6, 0.1, 19.6, 21.1, 19, 20, 20, 21.3, 21.6, 21.6, 0.2, 20.5, 0.25, 3.4)
a1c12 <- c(8.08, 8.08, 8.08, 8.2, 7.82, 8.7, 8.15, 7.4, 8.75, 7.89, 9.7, 7.6, 8.25, 7.57, 7.57, 7.57, 8.16, 8.16, 8.16, 7.29, 7, 7.25, 6.8)
sex12 <- c(45, 47, 46, 80, NA, 65, 32, 47, 66, 38, 72, 32, 72.5, 45.7, 50, 34.4, 47.5, 48.2, 48.9, 56.5, 46, 65.3, 50)





year12 <- scale(year12, scale = FALSE)
locna12 <- scale(locna12, scale = FALSE)
loceu12 <- scale(loceu12, scale = FALSE)
locas12 <- scale(locas12, scale = FALSE)
week12 <- scale(week12, scale = FALSE)
exena12 <- scale(exena12, scale = FALSE)
liradose12 <- scale(albi12, scale = FALSE)
rand12 <- scale(rand12, scale = FALSE) 
attr12 <- scale(attr12, scale = FALSE) 
detec12 <- scale(detec12, scale = FALSE)
repoth12 <- scale(repoth12, scale = FALSE)
dur12 <- scale(dur12, scale = FALSE)
a1c12 <- scale(a1c12, scale = FALSE)


c1o2d= data.frame(TE= c(c1o2$TE), seTE= c(c1o2$seTE), 
                  Publication_Year= year12,
                  Continent_NA_MN= locna12,
                  Continent_EU_MN= loceu12,
                  Continent_AS_MN= locas12,
                  Study_Duration= week12,
                  Med_Exena_Lira= exena12,
                  Liraglutide_Dose= liradose12,
                  Randomization_Bias= rand12,
                  Attrition_Bias= attr12,
                  Detection_Bias= detec12,
                  Reporting_Other_Bias= repoth12,
                  DM_Duration= dur12,
                  Baseline_A1c= a1c12)

# 
# c1o2.mi.test <- multimodel.inference(TE= "TE", seTE = "seTE", data = c1o2d, method='REML', test='knha',
#                                      predictors = c("Med_Exena_Lira", 
#                                                     "Med_Albi_Lira", 
#                                                     "Liraglutide_Dose",
#                                                     "Randomization_Bias",
#                                                     "Attrition_Bias",
#                                                     "Male_Sex", "DM_Duration", "Baseline_A1c"),
#                                      interaction=FALSE)
# c1o2.mi.test


c1o2.mi <- multimodel.inference(TE= "TE", seTE = "seTE", data = c1o2d, method='REML', test='knha',
                                predictors = c("Publication_Year", 
                                               "Continent_NA_MN", 
                                               "Continent_EU_MN", 
                                               "Continent_AS_MN", 
                                               "Study_Duration",
                                               "Med_Exena_Lira", 
                                               "Liraglutide_Dose",
                                               "Randomization_Bias",
                                               "Attrition_Bias",
                                               "Detection_Bias",
                                               "Reporting_Other_Bias",
                                               "DM_Duration", "Baseline_A1c"),
                                interaction=FALSE)
c1o2.mi



c1o2.mreg <- rma(yi = TE,
                 sei = seTE,
                 data = c1o2,
                 method = "REML",
                 mods = ~liradose12+a1c12+exena12+sex12,
                 test = "knha")
c1o2.mreg

c1o2.mreg1 <- rma(yi = TE,
                 sei = seTE,
                 data = c1o2,
                 method = "REML",
                 mods = ~liradose12+a1c12+exena12+sex12,
                 test = "knha")
c1o2.mreg1

anova(c1o2.mreg, c1o2.mreg1)



c1o3 <- metacr(GLPT1DMMETA_r,comp.no=1, outcome.no=3)
year13 <- c(2016, 2016, 2016, 2019, 2016, 2019, 2020, 2018, 2015, 2020, 2012, 2013, 2020, 2020, 2016, 2016, 2016, 2016, 2016, 2016, 2020, 2014, 2022)
#Location:   MN = 0, NA, Europe, Asia = 1
locna13 <- c(0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 1, 1, 1, 0, 0, 0, 0, 1, 0)
loceu13 <- c(0, 0, 0, 1, 1, 1, 1, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 1)
locas13 <- c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0)

week13 <- c(26, 26, 26, 26, 24, 52, 26, 24, 12, 26, 52, 52, 24, 26, 12, 12, 12, 26, 26, 26, 52, 52, 12)


#Regimen 0= lira 1= exena, albi
exena13 <- c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0)
albi13 <- c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0)
#Liradose 0= exena, albi  0.6, 1.2, 1.8= lira
liradose13 <- c(0.6, 1.2, 1.8, 1.8, 1.8, 1.8, 1.8, 1.8, 1.2, 1.8, 0.9, 0, 0, 0, 0.6, 1.2, 1.8, 0.6, 1.2, 1.8, 0, 0, 1.2)

#biases:   1= low risk, 2= medium, 3= high risk
rand13 <- c(1, 1, 1, 1, 1, 2, 1, 1, 1, 1, 3, 3, 2, 1, 2, 2, 2, 1, 1, 1, 1, 3, 2)
attr13 <- c(2, 2, 2, 3, 2, 2, 1, 1, 3, 3, 2, 1, 2, 3, 3, 3, 3, 2, 2, 2, 1, 3, 1)
detec13 <- c(1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2)
#Design:  Parallel =1 , Crossover = 2
des13 <- c(1, 1, 1, 1,  1, 1, 1, 2, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2)
repoth13 <- c(1, 1, 1, 1, 1, 2, 1, 2, 1, 1, 2, 2, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2)

dur13 <- c(20.9, 20.9, 21.1, 32.4, 22.5, 0.1, 20.5, 20.3, 18.9, 18.6, NA, 0.1, 19.6, 21.1, 19, 20, 20, 21.3, 21.6, 21.6, 0.2, 20.5, 3.4)
a1c13 <- c(8.08, 8.08, 8.08, 8.2, 8.7, 8.4, 8.15, 7.4, 8.75, 7.89, 7.65, 9.7, 7.6, 8.25, 7.57, 7.57, 7.57, 8.16, 8.16, 8.16, 7.29, 7, 6.8)
sex13 <- c(45, 47, 46, 80, 65, NA, 32, 47, 66, 38, NA, 72, 32, 72.5, 45.7, 50, 34.4, 47.5, 48.2, 48.9, 56.5, 46, 50)



year13 <- scale(year13, scale = FALSE)
locna13 <- scale(locna13, scale = FALSE)
loceu13 <- scale(loceu13, scale = FALSE)
locas13 <- scale(locas13, scale = FALSE)
week13 <- scale(week13, scale = FALSE)
exena13 <- scale(exena13, scale = FALSE)
liradose13 <- scale(albi13, scale = FALSE)
rand13 <- scale(rand13, scale = FALSE) 
attr13 <- scale(attr13, scale = FALSE) 
detec13 <- scale(detec13, scale = FALSE)
repoth13 <- scale(repoth13, scale = FALSE)
dur13 <- scale(dur13, scale = FALSE)
a1c13 <- scale(a1c13, scale = FALSE)


c1o3d= data.frame(TE= c(c1o3$TE), seTE= c(c1o3$seTE), 
                  Publication_Year= year13,
                  Continent_NA_MN= locna13,
                  Continent_EU_MN= loceu13,
                  Continent_AS_MN= locas13,
                  Study_Duration= week13,
                  Med_Exena_Lira= exena13,
                  Liraglutide_Dose= liradose13,
                  Randomization_Bias= rand13,
                  Attrition_Bias= attr13,
                  Detection_Bias= detec13,
                  Reporting_Other_Bias= repoth13,
                  DM_Duration= dur13,
                  Baseline_A1c= a1c13)

c1o3.mi <- multimodel.inference(TE= "TE", seTE = "seTE", data = c1o3d,
                                predictors = c("Publication_Year", 
                                               "Continent_NA_MN", 
                                               "Continent_EU_MN", 
                                               "Continent_AS_MN", 
                                               "Study_Duration",
                                               "Med_Exena_Lira", 
                                               "Liraglutide_Dose",
                                               "Randomization_Bias",
                                               "Attrition_Bias",
                                               "Detection_Bias",
                                               "Reporting_Other_Bias",
                                               "DM_Duration", "Baseline_A1c"),
                                interaction=FALSE)
c1o3.mi


.#
c1o3.mreg <- rma(yi = TE,
                 sei = seTE,
                 data = c1o3,
                 method = "ML",
                 mods = ~regimen13+sex13,
                 test = "knha")
c1o3.mreg
c1o3.mreg1 <- rma(yi = TE,
                  sei = seTE,
                  data = c1o3,
                  method = "ML",
                  mods = ~regimen13,
                  test = "knha")
c1o3.mreg1
anova(c1o3.mreg, c1o3.mreg1)




c1o4 <- metacr(GLPT1DMMETA_r,comp.no=1, outcome.no=4)
year14 <- c(2016, 2016, 2016, 2019, 2016, 2020, 2018, 2015, 2020, 2013, 2020, 2020, 2016, 2016, 2016, 2016, 2016, 2016, 2020, 2014, 2021, 2022)
#Location:   MN = 0, NA, Europe, Asia = 1
locna14 <- c(0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 1, 1, 1, 0, 0, 0, 0, 1, 0, 0)
loceu14 <- c(0, 0, 0, 1, 1, 1, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 1)
locas14 <- c(0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0)

week14 <- c(26, 26, 26, 26, 24, 26, 24, 12, 26, 52, 24, 26, 12, 12, 12, 26, 26, 26, 52, 52, 54, 12)


#Regimen 0= lira 1= exena, albi
exena14 <- c(0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0)
albi14 <- c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0)
#Liradose 0= exena, albi  0.6, 1.2, 1.8= lira
liradose14 <- c(0.6, 1.2, 1.8, 1.8, 1.8, 1.8, 1.8, 1.2, 1.8, 0, 0, 0, 0.6, 1.2, 1.8, 0.6, 1.2, 1.8, 0, 0, 1.8, 1.2)


#biases:   1= low risk, 2= medium, 3= high risk
rand14 <- c(1, 1, 1, 1, 1, 1, 1, 1, 1, 3, 2, 1, 2, 2, 2, 1, 1, 1, 1, 3, 1, 2)
attr14 <- c(2, 2, 2, 3, 2, 1, 1, 3, 3, 1, 2, 3, 3, 3, 3, 2, 2, 2, 1, 3, 3, 1)
detec14 <- c(1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2)
repoth14 <- c(1, 1, 1, 1, 1, 1, 2, 1, 1, 2, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2)

dur14 <- c(20.9, 20.9, 21.1, 32.4, 22.5, 20.5, 20.3, 18.9, 18.6, 0.1, 19.6, 21.1, 19, 20, 20, 21.3, 21.6, 21.6, 0.2, 20.5, 0.25, 3.4)
a1c14 <- c(8.08, 8.08, 8.08, 8.2, 8.7, 8.15, 7.4, 8.75, 7.89, 9.7, 7.6, 8.25, 7.57, 7.57, 7.57, 8.16, 8.16, 8.16, 7.29, 7, 7.25, 6.8)
sex14 <- c(45, 47, 46, 80, 65, 32, 47, 66, 38, 72, 32, 72.5, 45.7, 50, 34.4, 47.5, 48.2, 48.9, 56.5, 46, 65.3, 50)



year14 <- scale(year14, scale = FALSE)
locna14 <- scale(locna14, scale = FALSE)
loceu14 <- scale(loceu14, scale = FALSE)
locas14 <- scale(locas14, scale = FALSE)
week14 <- scale(week14, scale = FALSE)
exena14 <- scale(exena14, scale = FALSE)
liradose14 <- scale(albi14, scale = FALSE)
rand14 <- scale(rand14, scale = FALSE) 
attr14 <- scale(attr14, scale = FALSE) 
detec14 <- scale(detec14, scale = FALSE)
repoth14 <- scale(repoth14, scale = FALSE)
dur14 <- scale(dur14, scale = FALSE)
a1c14 <- scale(a1c14, scale = FALSE)


c1o4d= data.frame(TE= c(c1o4$TE), seTE= c(c1o4$seTE), 
                  Publication_Year= year14,
                  Continent_NA_MN= locna14,
                  Continent_EU_MN= loceu14,
                  Continent_AS_MN= locas14,
                  Study_Duration= week14,
                  Med_Exena_Lira= exena14,
                  Liraglutide_Dose= liradose14,
                  Randomization_Bias= rand14,
                  Attrition_Bias= attr14,
                  Detection_Bias= detec14,
                  Reporting_Other_Bias= repoth14,
                  DM_Duration= dur14,
                  Baseline_A1c= a1c14)

c1o4.mi <- multimodel.inference(TE= "TE", seTE = "seTE", data = c1o4d,
                                predictors = c("Publication_Year", 
                                               "Continent_NA_MN", 
                                               "Continent_EU_MN", 
                                               "Continent_AS_MN", 
                                               "Study_Duration",
                                               "Med_Exena_Lira", 
                                               "Liraglutide_Dose",
                                               "Randomization_Bias",
                                               "Attrition_Bias",
                                               "Detection_Bias",
                                               "Reporting_Other_Bias",
                                               "DM_Duration", "Baseline_A1c"),
                                interaction=FALSE)
c1o4.mi

c1o4.mreg <- rma(yi = TE,
                 sei = seTE,
                 data = c1o4,
                 method = "REML",
                 mods = ~locas14+liradose14+albi14+attr14,
                 test = "knha")
c1o4.mreg
c1o4.mreg1 <- rma(yi = TE,
                  sei = seTE,
                  data = c1o4,
                  method = "REML",
                  mods = ~sex14,
                  test = "knha")
c1o4.mreg1
anova(c1o4.mreg, c1o4.mreg1)
